function [M,VAR] = generic_nominal_moments(inflation_process,term0,t,termt,prob_NxYz, Phi_NxYz_t, grid_Y)
% computes the following:
% E[term(0)+term(t)] and Var[term(0)+term(t)]
% where the terms F(t)=term(0), term(t) are all of the form
% F(t) = fun_Nxz(N(t),x(t),z(t)) + const + slope_Y*Y(t) + slope_X*X(t) + slope_eps*eps(t)

    if t<=0; error('Expecting t>0'); end

    [NN, Nx, NY, Nz] = size(prob_NxYz);

    mean_X = inflation_process.mu_pi/(1 - inflation_process.rho_pi);
    var_X = (1 + 2*inflation_process.rho_pi*inflation_process.theta_pi ...
        + inflation_process.theta_pi^2)*(inflation_process.sigma_pi^2)...
        /(1 - inflation_process.rho_pi^2);
    var_eps = inflation_process.sigma_pi^2;
    cov_Xe = var_eps;
    
    cov_X0_Xt = (inflation_process.rho_pi^t)*var_X ...
        + (inflation_process.rho_pi^(t-1))*inflation_process.theta_pi*cov_Xe;
    cov_eps0_Xt = (inflation_process.rho_pi^t)*cov_Xe ...
        + (inflation_process.rho_pi^(t-1))*inflation_process.theta_pi*var_eps;
    
    Y_NxYz = permute(repmat(grid_Y(:),[1 NN Nx Nz]),[2 3 1 4]);

    % term0
    term0_NxYz = permute(repmat(term0.fun_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        + term0.slope_Y*Y_NxYz;
    mean_term0_NxYz = sum(prob_NxYz.*term0_NxYz,'all');
    var_term0_NxYz = sum(prob_NxYz.*(term0_NxYz.^2),'all') - mean_term0_NxYz^2;
    
    mean_term0_Xe = term0.const+ term0.slope_X*mean_X;
    var_term0_Xe = (term0.slope_X^2)*var_X + (term0.slope_eps^2)*var_eps ...
        + 2*term0.slope_X*term0.slope_eps*cov_Xe;

    % termt
    termt_NxYz = permute(repmat(termt.fun_Nxz,[1 1 1 NY]),[1 2 4 3]) ...
        + termt.slope_Y*Y_NxYz;
    mean_termt_NxYz = sum(prob_NxYz.*termt_NxYz,'all');
    var_termt_NxYz = sum(prob_NxYz.*(termt_NxYz.^2),'all') - mean_termt_NxYz^2;
    
    mean_termt_Xe = termt.const+ termt.slope_X*mean_X;
    var_termt_Xe = (termt.slope_X^2)*var_X + (termt.slope_eps^2)*var_eps ...
        + 2*termt.slope_X*termt.slope_eps*cov_Xe;

    % covariances
    cov_termt_term0_NxYz = sum(prob_NxYz(:).*term0_NxYz(:).*(Phi_NxYz_t*termt_NxYz(:))) ...
        - mean_termt_NxYz*mean_term0_NxYz;
    cov_termt_term0_Xe = termt.slope_X*term0.slope_X*cov_X0_Xt ...
        + termt.slope_X*term0.slope_eps*cov_eps0_Xt;

    % Mean
    M = mean_term0_NxYz + mean_termt_NxYz + mean_term0_Xe + mean_termt_Xe;

    % Var
    VAR = var_term0_NxYz + var_termt_NxYz + 2*cov_termt_term0_NxYz ...
        + var_term0_Xe + var_termt_Xe + 2*cov_termt_term0_Xe;

end